{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "from sympy import *\n",
    "init_printing()\n",
    "from IPython.display import display"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beam Segment Equations\n",
    "\n",
    "Loading along a beam is not always mathematically continuous.  Concentrated forces and moments, as well as starting and ending points of distributed loads interupt the mathematical continuity of the beam.  For this reason, it is convenient to segment the beam between discontinuities to form a series of mathematically continuous segments. This allows for the direct calculation of maximum/minimum forces, rotations, and deflections in each segment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Distributed Load Equation\n",
    "\n",
    "The equation for a linear distributed load is expressed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "remove_input"
    ]
   },
   "outputs": [],
   "source": [
    "w, w1, w2, x, L = symbols('omega, omega_1, omega_2, x, L')\n",
    "w = w1+(w2-w1)/L*x\n",
    "display(Eq(symbols('\\\\omega(x)'), w))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where:\n",
    "\n",
    "$\\omega_1=$ magnitude of distributed load at starting point\n",
    "\n",
    "$\\omega_2=$ magnitude of distributed laod at ending point\n",
    "\n",
    "$x_1=$ distributed load starting point\n",
    "\n",
    "$x_2=$ distributed load ending point"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shear Force Equation\n",
    "The change in the shear force is equal to the area under (or integral of) the load diagram. If the shear force at the start of the segment, $V_1$, is known, the shear force at any point, $x$, relative to the start of the segment can be expressed as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = symbols('V_1')\n",
    "V = V1 + integrate(w, (x, 0, x))\n",
    "display(Eq(symbols('V(x)'), V1 + Integral(symbols('\\\\omega(x)'), (x, 0, x))))\n",
    "display(Eq(symbols('V(x)'), V1 + Integral(w, (x, 0, x))))\n",
    "display(Eq(symbols('V(x)'), V))\n",
    "print(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Moment Equation\n",
    "\n",
    "Similarly, the change in the moment diagram is equal to the area under the shear diagram. If the moment at the start of the segment, $M_1$, is known, the moment at any point, $x$, relative to the start of the segment can be expressed as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M1 = symbols('M_1')\n",
    "M = M1 - integrate(V, (x, 0, x))\n",
    "display(Eq(symbols('M(x)'), M1 - Integral(symbols('V(x)'), (x, 0, x))))\n",
    "display(Eq(symbols('M(x)'), M1 - Integral(V, (x, 0, x))))\n",
    "display(Eq(symbols('M(x)'), M))\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Slope Equation\n",
    "The slope of the elastic curve, $\\theta$, can be obtained by integrating the moment diagram and dividing by the flexural rigidity, $EI$.\n",
    "\n",
    "$\\theta_1=$ slope (or rotation) at the start of the segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta1, E, I = symbols('theta_1, E, I')\n",
    "theta = theta1 - 1/(E*I)*integrate(M, (x, 0, x))\n",
    "display(Eq(symbols('\\\\theta(x)'), theta1 + Integral(symbols('M(x)'), (x, 0, x))))\n",
    "display(Eq(symbols('\\\\theta(x)'), theta1 + Integral(M, (x, 0, x))/(E*I)))\n",
    "display(Eq(symbols('\\\\theta(x)'), theta))\n",
    "print(theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deflection Equation\n",
    "Integrating one more time gives us the equation for the deflection along the segment.\n",
    "\n",
    "$\\Delta_1=$ Deflection at the start of the segment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Delta1 = symbols('Delta_1')\n",
    "Delta = Delta1 + integrate(theta, (x, 0, x))\n",
    "display(Eq(symbols('\\\\Delta(x)'), Delta1 + Integral(symbols('\\\\theta(x)'), (x, 0, x))))\n",
    "display(Eq(symbols('\\\\Delta(x)'), Delta1 + Integral(theta, (x, 0, x))))\n",
    "display(Eq(symbols('\\\\Delta(x)'), Delta))\n",
    "print(Delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the Equations\n",
    "Each of the equations derived above requires knowledge of the initial conditions at the start of the segment ($V_1$, $M_1$, $\\theta_1$, or $\\Delta_1$).  An analysis will need to be performed to determine these values at the start of the beam.\n",
    "\n",
    "The initial conditions at the start of the beam are the same as those at the start of the first segment.  If the initial conditions at the start of the first segment of the beam are known, the conditions at the start of the second segment can be found by calculating the conditions at the end of the first segment (from the equations with $x=L$) and adding any concentrated loads that may be at the junction between segments.  Likewise, the initial conditions for the third segment can be derived from the second segment, and so on until the last segment.  Any concentrated loads at the end of the last segment (including reactions) will need to be added to get shear and moment diagrams to close.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating $\\theta_1$\n",
    "\n",
    "For finite element solutions, the rotation $\\theta_1$ at the start of the beam might be known. However, if end releases have been applied to the member (e.g. members with hinges at one or both ends), then either $\\theta_1$ or $\\theta_2$ or both will not be known. These terms will be lost when the stiffness matrix is condensed to simulate the end releases. In this situation the slope-deflection equation may be used to solve for $\\theta_1$ as follows (since $M_1$, $M_2$, $\\Delta_1$, $\\Delta_2$, $FEM_1$, & $FEM_2$ will be known). The basic slope-deflection equation is:\n",
    "\n",
    "$M_{nf}=\\frac{2EI}{L}(2\\theta_n+\\theta_f-3\\psi)+FEM_{nf}$\n",
    "\n",
    "Substituting in appropriate values for each end of the segment yields the following equations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M1, M2, E, I, L, theta1, theta2, Delta1, Delta2 = symbols('M_1, M_2, E, I, L, theta_1, theta_2, Delta_1, Delta_2')\n",
    "FEM1, FEM2 = symbols('FEM_1, FEM_2')\n",
    "M1 = 2*E*I/L*(2*theta1 + theta2 - 3*(Delta2-Delta1)/L) + FEM1\n",
    "M2 = 2*E*I/L*(2*theta2 + theta1 - 3*(Delta2-Delta1)/L) + FEM2\n",
    "display(Eq(symbols('M_1'), M1))\n",
    "display(Eq(symbols('M_2'), M2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the equation $M_2$ slope-deflection equation for $\\theta_2$ gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta2 = solve(M2 - symbols('M_2'), theta2)[0]\n",
    "display(Eq(symbols('\\\\theta_2'), expand(theta2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Substituting $\\theta_2$ back into the slope-deflection equation for $M_1$ gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "M1 = M1.subs('theta_2', theta2)\n",
    "display(Eq(symbols('M_1'), expand(M1)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving this new equation for $\\theta_1$ gives:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta1 = solve(M1 - symbols('M_1'), theta1)[0]\n",
    "display(Eq(symbols('\\\\theta_1'), expand(theta1)))\n",
    "print(expand(theta1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finding Max and Min Values\n",
    "\n",
    "The points of max/min shear can be found by setting the derivative of the equation for shear equal to zero. Ther derivative of the shear force is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(Eq(Derivative(symbols('V'), x), diff(V, x)))\n",
    "display(Eq(symbols('x'), solve(diff(V, x), symbols('x'))[0]))\n",
    "print(solve(diff(V, x), symbols('x'))[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The points of max/min moment can be found by setting the derivative of the equation for moment equal to zero. The derivative of the moment is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(Eq(Derivative(symbols('M'), x), diff(M, x)))\n",
    "display(Eq(symbols('x_1'), solve(diff(M, x), symbols('x'))[0]))\n",
    "display(Eq(symbols('x_2'), solve(diff(M, x), symbols('x'))[1]))\n",
    "print(solve(diff(M, x), symbols('x'))[0])\n",
    "print(solve(diff(M, x), symbols('x'))[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
